/*
 *  thermoelastic.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 26/10/12.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "thermoelastic.h"

Thermoelastic::Thermoelastic()
:	Material(),
	matHyper(),
    T0(),
    a(),
    c0(),
    K0(),
    wk(),
    gam(),
    beta(),
    dsT(),
    G(),
	stVar()
{}

// ***

Thermoelastic::Thermoelastic(ifstream &in, const string& na, const string & t, Material* ma)
:	Material(na,t),
	matHyper(ma),
    T0(300),
    a(22.333E-5),
    c0(1830),
    K0(0.15),
    wk(0.0004),
    gam(2),
    beta(0.0),
    dsT(0.0),
    G(0.0),
    stVar(true)
{
	string keyword, line;
	vector<string>	conductivityParam;
	stringstream stnum;
	
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("EXPASION"))
        {
			in >> a;
        }
        if(!keyword.compare("REFTEMP"))
        {
        	in >> T0;
        }
        if(!keyword.compare("SPECIFIC-HEAT"))
        {
        	in >> c0;
        }
		if(!keyword.compare("CONDUCTIVITY"))
		{
			getline(in,line);
			istringstream iss(line);
			while (iss >> keyword )
			{
				//cout << "keyword is " << keyword << endl;
				conductivityParam.push_back(keyword);
			}
			stnum << conductivityParam[0];
			stnum >> K0;

			FourierType = conductivityParam[1];
		}
		if (!keyword.compare("SOFTENING"))
		{
			in >> wk;
		}
		if (!keyword.compare("BULKPARAM"))
        {
            in >> gam;
        }
		if (!keyword.compare("LINEAREXPANSION"))
		{
		    in >> beta;
		}
		if (!keyword.compare("ENTROPY"))
        {
            stVar = false;
        }
    }

    hypertype = matHyper->hypertype;
    density = matHyper->getDensity();

    Ci.setZero(6);
    dG.setZero(6);
    dcT.setZero(6);
}       

// ***

void Thermoelastic::calcW( const VectorXd & C , const double & T )
{
	matHyper->calcW(C); // Here the invariants are computed for C
	calcG(matHyper->getJ());

	if(!hypertype.compare("coupled"))
    {
		W = matHyper->getW() * T / T0 - 3.0 * a * G * (T - T0) + c0 * (T - T0 - T * log(T / T0));
    }
    else
    {
    	W = matHyper->getWiso() * T / T0 - 3.0 * a * K * G * (T - T0) + c0 * (T - T0 - T * log(T / T0));
    }
}

// ***

void Thermoelastic::calcU(const VectorXd & C, const double & ent)
{
	calcTemp(C,ent);
	calcW(C,Temp);
        
	U = W +  Temp * ent;
}

// ***

void Thermoelastic::calcEntropy(const VectorXd & C, const double & T)
{
	matHyper->calcW(C); // Here the invariants are computed for C
	calcG(matHyper->getJ());

	if(!hypertype.compare("coupled"))
	{
		Ent = - (matHyper->getW() / T0 - 3.0 * a * G - c0 * log(T / T0));
	}
	else
	{
		Ent = - (matHyper->getWiso() / T0 - 3.0 * a * K * G - c0 * log(T / T0));
	}
}

// ***

void Thermoelastic::calcTemp(const VectorXd & C, const double & ent)
{
	matHyper->calcW(C); // Here the invariants are computed for C

	if(!hypertype.compare("coupled"))
    {
		Temp = T0 * exp((ent + matHyper->getW() / T0 - 3.0 * a * K * log(matHyper->getJ()) ) / c0);
    }
    else
    {
    	Temp = T0 * exp((ent + matHyper->getWiso() / T0 - 3.0 * a * K * log(matHyper->getJ()) ) / c0);
    }
}

// ***

void Thermoelastic::calcS( const VectorXd & C, const double & T )
{
	matHyper->calcS(C);

	switch(stVar)
	{
		case true:

			if(!hypertype.compare("coupled"))
			{
				S = T / T0 * matHyper->getS() - 3.0 * a * matHyper->getS();
			}
			else
			{
				calcDG(C);
				S = T / T0 * matHyper->getSiso() - 3.0 * a * K * (T - T0) * dG;
			}
			break;

		case false:
			// In this case T = entropy
			calcTemp(C,T);
			calcW(C,Temp);


			if(!hypertype.compare("coupled"))
			{
				calcDcTemp(C,T);
				S = Temp / T0 * matHyper->getS() + getW() / T0 * dcT - 3.0 * a * matHyper->getS();
			}
			else
			{
				calcDG(C);
				S = Temp / T0 * matHyper->getSiso() - 3.0 * a * K * (Temp - T0) * dG;
			}
			break;
	}
}

// ***

void Thermoelastic::calcCt( const VectorXd & C, const double & T)
{

}

// ***

void Thermoelastic::calcDcTemp(const VectorXd & C, const double & ent)
{
	calcTemp(C,ent);
	calcS(C,T0);
	calcDG(C);

	if(!hypertype.compare("coupled"))
    {
		dcT = Temp / c0 * ( S / 2.0 / T0 - 3.0 * a * K0 * dG );
    }
    else
    {
    	dcT = Temp / c0 * ( Siso / 2.0 / T0 - 3.0 * a * K0 * dG );
    }
}

// ***

void Thermoelastic::calcDsTemp(const VectorXd & C, const double & ent)
{
	calcTemp(C,ent);
    dsT = Temp / c0;
}

// ***


void Thermoelastic::calcG(const double & J )
{
	if (beta == 0.0)
	{
		G = 1.0 / gam * pow((J - 1),gam);
	}
	else
	{
		//To be corrected
		G = - 3.0 * beta * matHyper->getW();
	}
}

// ***

void Thermoelastic::calcDG(const VectorXd & C )
{
	matHyper->calcInvariants(C);
	matHyper->calcCiTensor(C);
	VoigtFormat(matHyper->getCiTensor(), Ci);
        
    dG = pow((matHyper->getJ() - 1), gam - 1.0) * matHyper->getJ() / 2 * Ci ;
}

// ***

void Thermoelastic::calcConductivity(const double & T)
{
	// Conductivity
	if(!FourierType.compare(""))
	{
		K = K0 * (1.0 - wk * (T - T0));
	}
	else if (!FourierType.compare("log"))
	{
		K = - K0 / T * (1.0 - wk * (T - T0));
	}
	else if (!FourierType.compare("inv"))
	{
		K = K0 * pow(T,2.0) * (1.0 - wk * (T - T0));
	}
}

// ***

double Thermoelastic::getT0()
{
	// Reference Temperature
	return T0;
}

// ***

double Thermoelastic::geta()
{
	// Linear expansion coefficient
	return a;
}

// ***

double Thermoelastic::getc0()
{
	// Specific heat capacity
	return c0;
}

// ***

double Thermoelastic::getK0()
{
	// Conductivity at T = T0
	return K0;
}

// ***

double Thermoelastic::getK()
{
	// Conductivity at T = T0
	return K;
}

// ***

double Thermoelastic::getwk()
{
	// Softening parameter
	return wk;
}

// ***

bool Thermoelastic::getMode()
{
	// State Variable: T or s?
	return stVar;
}
// ***
